%   Bestim.m
% - results calculation by cells is implemented in file "Bestimcell.m"
%
% %%%%%% Source of results %%%%%%%%%%
% - 
% - sample sizes from 5_sample2_cells_stats.xlsx
% - transition parameters from 5_sample2_wholesampleAndcells_estimates_Censoring
% - wage distribution from 5_sample2_wages_distribution.xls
%          which summarizes different coal vs. non-coal parameters across
%          groups
% 

clear; clc;

%% Baseline or robustness checks
 sample=2 % baseline
% sample=7 % using only post-2000 data
% sample=10 % robustness test for RR: doubling separation rates in NC sectors
% the other samples (8,9) cut data into different cells - use Bestimcell

%% estimation results %%
%% Monthly transition parameters & moments of monthly wages

if (sample==2 || sample==10)
%% Baseline values %%
% updated 16th Dece 2022 with run from 12th Dec
% - transition parameters from 5_sample2_wholesampleAndcells_estimates_Censoring
base.deltaC=30.5*0.000075719710;
    fprintf('every X months people are fired in coal = %f\n',  1/base.deltaC)
base.deltaNC=30.5*0.000440721354;
    fprintf('every X months people are fired in non-coal = %f\n', 1/base.deltaNC)
    if sample==10
        base.deltaNC=base.deltaNC*2;
    end
base.lamC=30.5*0.000158231643;
    fprintf('every X months unemp coal-workers find coal job = %f\n', 1/base.lamC)
base.lamNC=30.5*0.000909964725; 
    fprintf('every X months unemp coal-workers find non-coal job = %f\n', 1/base.lamNC)
base.lamZC=30.5*0.001812109028;
    fprintf('every X months unemp non-coal-workers find job = %f\n', 1/base.lamZC)
base.rho=30.5*0.000069348071;
    fprintf('every X months employed people retire = %f\n', 1/base.rho)
    fprintf('every X years employed people retire = %f\n', (1/base.rho)/12)
% - wage distribution from 5_sample2_wages_distribution.xls
base.wCmu  = 7.652580513055; %wcoal_lognormal
    fprintf('exp(mu)=median of lognormally distributed wages in coal job = %f\n', exp(base.wCmu))
base.wNCmu= 7.325744710244; 
    fprintf('exp(mu)=median of lognormally distributed wages in non-coal job= %f\n', exp(base.wNCmu))
base.wCsig =0.307794496258; 
base.wNCsig = 0.325125904498;   
elseif sample==7
    base.rho   = 30.5.*[0.000049773599];
    base.deltaC= 30.5.*[0.000037578119];
    base.deltaNC=30.5.*[0.000418117322];
    base.lamC  = 30.5.*[0.000215728440];
    base.lamNC = 30.5.*[0.001132513092];
    base.lamZC = 30.5.*[0.001901880197];
    base.wCmu  = 7.652575171370; 
    base.wNCmu= 7.324437497894; 
    base.wCsig = 0.505323829299 ; 
    base.wNCsig = 0.384629126276;   
% We are using 1st percentile as min and 99th percentile as max
% minfrac: since using only for integral - no harm in going assuming
% possible lower/higher wage offers (but maybe avoid extreme values of actual
% min/max)
end
% 
%% calculation of standard errors of transition parameters (whole pop)
%
% all these from 
% TransitionParametersMainResults.xlsx (below in sheet ByCell - Samples) 
%     &  5_sample2_wages_distribution.xlsx (sheet "wcoal_normal")
%
if (sample==2 ||sample==10)
    T.rho = 24148;
    T.lamC = 3876;
    T.lamNC = 22275;
    T.deltaC = T.lamC + T.lamNC;
    T.lamZC =  171764;
    T.deltaNC = 176256;
    T.wC = 22272 ;      % taken from 5_sample2_wages_distribution.xlsx
    T.wNC = 15391 ;     % taken from 5_sample2_wages_distribution.xlsx
elseif sample==7
    T.rho =     4970;
    T.lamC =    881;
    T.lamNC =   4625;
    T.deltaC =  3688; 
    T.lamZC =   118606;
    T.deltaNC=  109886;
    % sample sizes below taken from 5_sample7_wages_distribution (sum_w)
    T.wC =      2863;
    T.wNC=      2930;        
end

[error.rho,     cilo.rho,   cihi.rho]       = standarderror(base.rho,T.rho);
[error.deltaC,  cilo.deltaC,cihi.deltaC]    = standarderror(base.deltaC,T.deltaC);
[error.deltaNC, cilo.deltaNC,cihi.deltaNC]  = standarderror(base.deltaNC,T.deltaNC);
[error.lamC,    cilo.lamC,  cihi.lamC]      = standarderror(base.lamC,T.lamC);
[error.lamNC,   cilo.lamNC, cihi.lamNC]     = standarderror(base.lamNC,T.lamNC);
[error.lamZC,   cilo.lamZC, cihi.lamZC]     = standarderror(base.lamZC,T.lamZC);
error.wC= base.wCsig./sqrt(T.wC);
error.wNC=base.wNCsig./sqrt(T.wNC);

rlabse = {'s.e.rho', 'CI-high rho', 'CI-low rho','s.e.deltaC', 'CI-hi deltaC', 'CI-lo deltaC','s.e.deltaNC','CI-hi deltaNC','CI-lo deltaNC', 's.e.lamC', 'CI-hi lamC', 'CI-lo lamC', 's.e.lamNC', 'CI-hi lamNC', 'CI-lo lamNC','s.e.lamZC', 'CI-hi lamZC', 'CI-lo lamZC','s.e.wC', 's.e.wNC'};
clabse= {'Baseline model (whole pop)'};
    matse(1,1) = error.rho;
    matse(2,1) = cilo.rho;
    matse(3,1) = cihi.rho;
    matse(4,1) = error.deltaC;
    matse(5,1) = cilo.deltaC;
    matse(6,1) = cihi.deltaC;
    matse(7,1) = error.deltaNC;
    matse(8,1) = cilo.deltaNC;
    matse(9,1) = cihi.deltaNC;
    matse(10,1) = error.lamC;
    matse(11,1) = cilo.lamC;
    matse(12,1) = cihi.lamC;
    matse(13,1) = error.lamNC;
    matse(14,1) = cilo.lamNC;
    matse(15,1) = cihi.lamNC;
    matse(16,1) = error.lamZC;
    matse(17,1) = cilo.lamZC;
    matse(18,1) = cihi.lamZC;
    matse(19,1) = error.wC;
    matse(20,1) = error.wNC;

% output results to latex
path= '../../../paper/' ;
if sample==2
    file='transebase.tex';
elseif sample==7
    file='transebas7.tex';
elseif sample==10
    file='transebas10.tex';
end
filename= [path file];
matrix2latex(matse, filename, 'rowLabels', rlabse, 'columnLabels', clabse, 'alignment', 'c', 'format', '%-6.8f', 'size', 'small');
        
    %% Different groups for welfare simulation: baseline vs. counterfactuals 

    % group1: baseline results for welfare costs of coal job loss (WFC)
    % group2: counterf. WFC with F(w)^NC = F(w)^C
    % group3: counterf. WFC with delta^NC = delta^C
    % group4: counterf. WFC with (i) F(w)^NC = F(w)^C; (ii) delta^NC = delta^C
Nx=4;

% - using monthly transition rates from daily transition rates
param.deltaC= base.deltaC*ones(1,Nx);
param.deltaNC=base.deltaNC*ones(1,Nx);
param.lamC =  base.lamC*ones(1,Nx);
param.lamNC=  base.lamNC*ones(1,Nx);
param.rho  =  base.rho*ones(1,Nx); 
param.wCmu  = base.wCmu*ones(1,Nx) ;  
param.wNCmu = base.wNCmu*ones(1,Nx) ; 
param.wCsig = base.wCsig*ones(1,Nx);  
param.wNCsig= base.wNCsig*ones(1,Nx); 
%{
param.wCmin = base.wCmin*ones(1,Nx) ;  
param.wCmax = base.wCmax*ones(1,Nx) ;  
param.wNCmin = base.wNCmax*ones(1,Nx) ;  
param.wNCmax = base.wNCmax*ones(1,Nx) ;  
%}
% Group1 is for the baseline values 

% Group 2 tests WFC if NC wages as large as C (same w-distribution)
param.wNCmu(2)  = base.wCmu; 
param.wNCsig(2) = base.wCsig; 
%{
param.wNCmin(2) = base.wCmin; 
param.wNCmax(2) = base.wCmax; 
%}
% Group 3 tests WFC if job-loss rate in NC = C (same deltas)
param.deltaNC(3) = base.deltaC;

% Group 4 tests WFC with F(w^C)=F(w^NC) (group2) & same deltas (group3)
param.deltaNC(4)= base.deltaC;
param.wNCmu(4)  = base.wCmu; 
param.wNCsig(4) = base.wCsig; 
%{
param.wNCmin(4) = base.wCmin; 
param.wNCmax(4) = base.wCmax; 
%}

%% Now for calculation of values of unemployment, coal job & welfare cost:
%% call Calgo.m

for fut=0:1
    if fut==1 
        if sample==2
        data.label= 'resbase.tex'; % label latex file with results produced by Calgo
        elseif sample==7
        data.label= 'resbas7.tex'; % label latex file with results produced by Calgo
        elseif sample==10
        data.label= 'resbas10.tex'; % label latex file with results produced by Calgo
        end
    Calgo
    elseif fut==0
    data.label= 'resplusc.tex'; % label latex file with results produced by Calgo
    Calgo
    end
end

%% More analysis of results

% Decomposing baseline welfare costs:
wfc_base=matres(1,3);
% (1) Fraction due to differences in wages in C vs. NC sectors
wfc_samew=matres(2,3);
frac1=(wfc_base-wfc_samew)/wfc_base;
fprintf('fraction of welfare costs due to differences in wages = %f\n', frac1)
% (2) Fraction due to differences in jo loss prob across C vs. NC sectors
wfc_samedelta=matres(3,3);
frac2=(wfc_base-wfc_samedelta)/wfc_base;
fprintf('fraction of welfare costs due to differences in job loss prob = %f\n', frac2)
% (3) Fraction due to both differences in wages & job loss probs 
wfc_samewanddelta=matres(4,3);
frac3=(wfc_base-wfc_samewanddelta)/wfc_base;
fprintf('fraction of welfare costs due to diffs in wage & job loss prob = %f\n', frac3)


